suppressPackageStartupMessages({
library(vegan)
library(here)
library(magrittr)
library(cowplot)
library(sf)
library(tmap)
library(tidyverse)
})
kelp <- read.csv(here("data", "Peces_KelpForest_2011-2013.csv"),
stringsAsFactors = F)
kelp11 <- kelp %>%
filter(year == 2011)
Para datos 2011
kelp11 %>%
group_by(location, site, level, transect) %>%
tally() %>%
ungroup() %>%
select(-n) %>%
group_by(location, site, level) %>%
tally() %>%
knitr::kable(caption = "Numero de transectos por localidad, sitio (N-S) y nivel (F-MA), incluyendo registros de 'fuera de transecto'.")
| location | site | level | n |
|---|---|---|---|
| ASA | N | Bottom | 4 |
| ASA | N | Midwater | 4 |
| ASA | S | Bottom | 4 |
| ASA | S | Midwater | 3 |
| BMA | N | Bottom | 4 |
| BMA | N | Midwater | 4 |
| BMA | S | Bottom | 4 |
| BMA | S | Midwater | 4 |
| CKE | N | Bottom | 9 |
| CLO | N | Bottom | 3 |
| CLO | N | Midwater | 2 |
| CLO | S | Bottom | 4 |
| CLO | S | Midwater | 4 |
| COL | N | Bottom | 1 |
| COL | S | Bottom | 2 |
| ERE | N | Bottom | 4 |
| ERE | N | Midwater | 4 |
| ERE | S | Bottom | 4 |
| ERE | S | Midwater | 2 |
| ERO | N | Bottom | 4 |
| ERO | N | Midwater | 4 |
| ERO | S | Bottom | 4 |
| ERO | S | Midwater | 4 |
| ISME | N | Bottom | 4 |
| ISME | N | Midwater | 4 |
| ISME | S | Bottom | 4 |
| ISME | S | Midwater | 4 |
| ITSP | N | Bottom | 4 |
| ITSP | N | Midwater | 4 |
| ITSP | S | Bottom | 4 |
| ITSP | S | Midwater | 4 |
| PCH | S | Bottom | 4 |
| PCH | S | Midwater | 4 |
| RET | N | Bottom | 4 |
| RET | N | Midwater | 4 |
| RET | S | Bottom | 4 |
| RET | S | Midwater | 3 |
| SMI | N | Bottom | 4 |
| SMI | N | Midwater | 4 |
| SMI | S | Bottom | 4 |
| SMI | S | Midwater | 4 |
| SSI | N | Bottom | 4 |
| SSI | N | Midwater | 4 |
| SSI | S | Bottom | 4 |
| SSI | S | Midwater | 4 |
| STO | N | Bottom | 3 |
| STO | N | Midwater | 3 |
| STO | S | Bottom | 3 |
| STO | S | Midwater | 3 |
| VTR | N | Bottom | 3 |
| VTR | N | Midwater | 2 |
| VTR | S | Bottom | 1 |
| VTR | S | Midwater | 1 |
De la tabla anterior vemos que Campo Kenedy, Colonet, Campo Lopez, Eréndira y Vaye Tranquilo no tienen suficientes muestras (mínimo 3 transectos por nivel por zona), por lo que las excluimos en el siguiente análisis.
Todas las matrices de distancias son calculadas por distancia de Bray-Curtis.
Necesitamos una matriz de presencia / ausencia por especie para correr el ANOSIM entre F y MA
localidades_excluidas <- c("CLO", "COL", "ERE", "VTR", "CKE", "PCH")
data_test1 <- kelp11 %>%
filter(!location %in% localidades_excluidas) %>%
group_by(location, site, level, latitude, longitude, genusspecies) %>%
summarize(abundance = sum(abundance, na.rm = T)) %>%
ungroup() %>%
mutate(abundance = 1) %>%
spread(genusspecies, abundance, fill = 0)
data_test1_groups <- data_test1 %>%
select(location, site, level) %>%
mutate(loc_site = paste(location, site, sep = "-"))
data_test1_samples <- data_test1 %>%
select(-c(location, site, level, latitude, longitude)) %>%
vegdist(method = "bray")
Habiendo filtrado las localidades que no tenian suficientes muestras, obtenemos la siguiente lista de localidades: ASA, BMA, ERO, ISME, ITSP, RET, SMI, SSI, STO.
Los datos de arriba estan listos para correr los ANOSIMs. El objeto data_test1_groups tiene los factores o tratamientos correspondientes a cada combinacion localidad - sitio - nivel. Los datos tienen la suma a traves de los 4 transectos de cada localidad - sitio - nivel, y son convertidos a presencia / ausencia.
A continuación los resultados de ADONIS (ANOSIMS con ANOVAS de matrices de distancia), que son más robustos para anidaciones (diferente a dos vías).
Reconociendo que las diferencias entre Fondo y MA pueden ser generadas por procesos especificos a cada localidad, corremos el ANOSIM anidando los transectos F y MA dentro de localidad - sitio. En este caso, el ANOSIM anidado indica que dentro de cada localidad - sitio, existen diferencias entre el fondo y media agua.
perm <- how(nperm = 999)
setBlocks(perm) <- data_test1_groups$loc_site
set.seed(43)
adonis(formula = data_test1_samples ~ level, data = data_test1_groups, permutations = perm)
##
## Call:
## adonis(formula = data_test1_samples ~ level, data = data_test1_groups, permutations = perm)
##
## Blocks: data_test1_groups$loc_site
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## level 1 0.1391 0.13906 1.2869 0.03647 0.01 **
## Residuals 34 3.6740 0.10806 0.96353
## Total 35 3.8131 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
La comparacion similar anidando a nivel de localidad (ignorando las diferencias entre N y S) indica que dentro de cada localidad, independiente si es sitio N o S, también hay diferencias entre F y MA.
setBlocks(perm) <- data_test1_groups$location
set.seed(43)
(test1a <- adonis(formula = data_test1_samples ~ level, data = data_test1_groups, permutations = perm))
##
## Call:
## adonis(formula = data_test1_samples ~ level, data = data_test1_groups, permutations = perm)
##
## Blocks: data_test1_groups$location
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## level 1 0.1391 0.13906 1.2869 0.03647 0.034 *
## Residuals 34 3.6740 0.10806 0.96353
## Total 35 3.8131 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sin embargo, podemos correr tambien la comparacion de todos los fondos contra todos los MA, ignorando la variabilidad explicada por cada localidad - sitio. En este caso, las comunidades (determinadas por presencia / ausencia) no son suficientemente diferentes.
set.seed(43)
adonis(formula = data_test1_samples ~ level, data = data_test1_groups)
##
## Call:
## adonis(formula = data_test1_samples ~ level, data = data_test1_groups)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## level 1 0.1391 0.13906 1.2869 0.03647 0.268
## Residuals 34 3.6740 0.10806 0.96353
## Total 35 3.8131 1.00000
Ahora debenos comparar entre sitios, anidando a nivel de localidad. Es decir, comparar dentro de cada localidad los sitios N contra S, sin tomar en cuenta la variación que hay entre F y MA. Los resultados indican que para cada localidad, N y S (juntando F y MA) son suficientemente similares como para considerar cada localidad como un grupo homogéneo. No hacemos la comparación N - S sin anidar, pues N - S solamente se definen así dentro de cada localidad, pero no representan fuentes sistemáticas de variación dentro de cada sitio. La historia sería diferente si fuera Expuesto - Protegido, pero eso es en los datos de 2013, no en estos.
setBlocks(perm) <- data_test1_groups$location
set.seed(43)
adonis(formula = data_test1_samples ~ site, data = data_test1_groups, permutations = perm)
##
## Call:
## adonis(formula = data_test1_samples ~ site, data = data_test1_groups, permutations = perm)
##
## Blocks: data_test1_groups$location
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## site 1 0.0968 0.09684 0.88599 0.0254 0.195
## Residuals 34 3.7162 0.10930 0.9746
## Total 35 3.8131 1.0000
Finalmente, comparamos entre localidades. En este caso no hay anidacion, pues la localidad es el elemento jerárquico más alto. Los resultados de ADONIS indican diferencias significativas entre sitios.
set.seed(43)
adonis(formula = data_test1_samples ~ location, data = data_test1_groups)
##
## Call:
## adonis(formula = data_test1_samples ~ location, data = data_test1_groups)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## location 8 1.9286 0.241076 3.4541 0.50579 0.001 ***
## Residuals 27 1.8845 0.069795 0.49421
## Total 35 3.8131 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Necesitamos obtener la densidad promedio de cada especie a nivel de localidad - sitio - nivel. En este caso es necesario filtrar los datos de “fuera de transecto”, o transecto 0. Las densidades promedio son transformadas por raiz cuadrada.
data_test2 <- kelp11 %>%
filter(transect > 0,
!location %in% localidades_excluidas) %>%
group_by(location, site, level, transect, genusspecies) %>%
summarize(abundance = sum(abundance)) %>%
ungroup() %>%
group_by(location, site, level, genusspecies) %>%
summarize(abundance = mean(abundance, na.rm = T)) %>%
ungroup() %>%
mutate(abundance = sqrt(abundance)) %>%
spread(genusspecies, abundance, fill = 0)
data_test2_groups <- data_test2 %>%
select(location, site, level) %>%
mutate(loc_site = paste(location, site, sep = "-"))
data_test2_samples <- data_test2 %>%
select(-c(location, site, level)) %>%
vegdist(method = "bray")
set.seed(43)
mds <- metaMDS(data_test2_samples, trace = F)
stress <- paste("2D Stress =", formatC(mds$grstress, digits = 4, format = "f"))
p1 <- cbind(data_test2_groups, scores(mds)) %>%
ggplot(aes(x = NMDS1, y = NMDS2, sitio = site)) +
geom_point(size = 4, aes(color = location, shape = level)) +
coord_equal() +
scale_color_brewer(palette = "Paired") +
annotate(geom = "text", x = 0, y = 0.5, label = stress)
p1
Version animada:
plotly::ggplotly(p1)
El nMDS muestra claramente diferencias en la densidades entre F y MA.
Podemos corroborar haciendo un nMDS con todos los transectos (sin promediar). En este caso, cada localidad deberia de tener al rededor de 12 puntos (3 NF, 3NM, 3SF, 3SM). Las siguientes figuras es exactamente la misma informacion, pero con diferente representacion. La primera solamente tiene informacion visual de nivel (F vs MA), la segund aincluye tambien color segun localidad.
data_test2 <- kelp11 %>%
filter(transect > 0,
!location %in% localidades_excluidas) %>%
group_by(location, site, level, transect, genusspecies) %>%
summarize(abundance = sum(abundance)) %>%
ungroup() %>%
mutate(abundance = sqrt(abundance)) %>%
spread(genusspecies, abundance, fill = 0)
data_test2_groups <- data_test2 %>%
select(location, site, level) %>%
mutate(loc_site = paste(location, site, sep = "-"))
data_test2_samples <- data_test2 %>%
select(-c(location, site, level)) %>%
vegdist(method = "bray")
set.seed(43)
mds <- metaMDS(data_test2_samples, trace = F)
stress <- paste("2D Stress =", formatC(mds$grstress, digits = 4, format = "f"))
p1 <- cbind(data_test2_groups, scores(mds)) %>%
ggplot(aes(x = NMDS1, y = NMDS2, sitio = site)) +
geom_point(size = 4, aes(shape = level, color = level)) +
coord_equal() +
scale_color_brewer(palette = "Paired") +
annotate(geom = "text", x = 0, y = 0.5, label = stress)
p2 <- cbind(data_test2_groups, scores(mds)) %>%
ggplot(aes(x = NMDS1, y = NMDS2, sitio = site)) +
geom_point(size = 4, aes(color = location, shape = level)) +
coord_equal() +
scale_color_brewer(palette = "Paired") +
annotate(geom = "text", x = 0, y = 0.5, label = stress)
plot_grid(p1, p2, ncol = 1, labels = "AUTO")
El ANOSIM confirma las diferencia sobservadas en el nMDS de arriba. Este anosim compara las diferencias entre fondo y media agua, anidando por localidad y sitio
perm <- how(nperm = 999)
setBlocks(perm) <- data_test2_groups$loc_site
set.seed(43)
adonis(formula = data_test2_samples ~ level, data = data_test2_groups, permutations = perm)
##
## Call:
## adonis(formula = data_test2_samples ~ level, data = data_test2_groups, permutations = perm)
##
## Blocks: data_test2_groups$loc_site
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## level 1 3.2863 3.2863 17.517 0.14656 0.001 ***
## Residuals 102 19.1356 0.1876 0.85344
## Total 103 22.4219 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Usando solamente los datos de FONDO, corremos un ANOSIM de dos vias comparando entre localidad y sitio. En este caso identificamos que no hay diferencias entre sitio, pero si entre localidades.
# Filtramos los datos para tener solamente FONDO
data_test2_groups_B <- data_test2 %>%
select(location, site, level) %>%
filter(level == "Bottom") %>%
mutate(loc_site = paste(location, site, sep = "-"))
data_test2_samples_B <- data_test2 %>%
filter(level == "Bottom") %>%
select(-c(location, site, level)) %>%
vegdist(method = "bray")
set.seed(43)
adonis(formula = data_test2_samples_B ~ site + location, data = data_test2_groups_B)
##
## Call:
## adonis(formula = data_test2_samples_B ~ site + location, data = data_test2_groups_B)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## site 1 0.2028 0.20279 1.7547 0.02111 0.086 .
## location 8 4.3178 0.53973 4.6702 0.44951 0.001 ***
## Residuals 44 5.0850 0.11557 0.52938
## Total 53 9.6056 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Usando solamente los datos de MEDIA AGUA, corremos un ANOSIM de dos vias comparando entre localidad y sitio. En este caso, las diferencias entre localidades son significativas, pero no las diferencias entre sitio.
# Filtramos los datos para tener solamente MEDIA AGUA
data_test2_groups_M <- data_test2 %>%
select(location, site, level) %>%
filter(level == "Midwater") %>%
mutate(loc_site = paste(location, site, sep = "-"))
data_test2_samples_M <- data_test2 %>%
filter(level == "Midwater") %>%
select(-c(location, site, level)) %>%
vegdist(method = "bray")
set.seed(43)
adonis(formula = data_test2_samples_M ~ site + location, data = data_test2_groups_M)
##
## Call:
## adonis(formula = data_test2_samples_M ~ site + location, data = data_test2_groups_M)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## site 1 0.2619 0.26189 1.7431 0.02748 0.100 .
## location 8 3.2584 0.40731 2.7110 0.34192 0.001 ***
## Residuals 40 6.0096 0.15024 0.63060
## Total 49 9.5300 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Ahora comparamos la composicion de la comunidad entre niveles y entre localidades. Esto nos indica que hay diferencias entre ambos.
set.seed(43)
adonis(formula = data_test2_samples ~ level + location, data = data_test2_groups)
##
## Call:
## adonis(formula = data_test2_samples ~ level + location, data = data_test2_groups)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## level 1 3.2863 3.2863 21.8719 0.14656 0.001 ***
## location 8 5.0121 0.6265 4.1698 0.22354 0.001 ***
## Residuals 94 14.1235 0.1502 0.62990
## Total 103 22.4219 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Contestamos ambas preguntas con el IVB y SIMPER
source('~/GitHub/bvi/bvi.R')
source('~/GitHub/bvi/bvi_plot.R')
source('~/GitHub/bvi/bvi_col.R')
source('~/GitHub/bvi/bvi_boxplot.R')
Para el IVB, usamos unicamente los datos del 2011 y con las mismas localidades que hemos usado hasta ahora. (Al incluir todas las localidades del 2011 los resultados son los mismos). Calculamos las densidades promedio por localidad y especie. Cuando una especie no estaba presente, recibe densidad de 0 (Esto se puede cambiar, pero los resultados son los mismos).
site_data <- kelp %>%
group_by(year, location, site, latitude, longitude) %>%
tally()
bvi_results <- filter(kelp, transect > 0) %>%
filter(year == 2011,
!location %in% localidades_excluidas) %>%
group_by(location, site, level, transect, genusspecies) %>%
summarize(n = sum(abundance)) %>%
spread(location, n, fill = 0) %>%
gather(location, n,-c(site, level, transect, genusspecies)) %>%
group_by(location, genusspecies) %>%
summarize(n = mean(n)) %>%
ungroup() %>%
select(location, genusspecies, n) %>%
spread(location, n) %>%
rename(Spp = genusspecies) %>%
bvi(sum = F, others = T)
knitr::kable(bvi_results)
| Spp | ASA | BMA | ERO | ISME | ITSP | RET | SMI | SSI | STO | BVI | rBVI |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Oxyjulis californica | 17 | 16 | 11 | 17 | 17 | 15 | 17 | 16 | 17 | 143 | 10.206995 |
| Semicossyphus pulcher | 16 | 15 | 12 | 16 | 16 | 14 | 14 | 14 | 15 | 132 | 9.421841 |
| Chromis punctipinnis | 12 | 17 | 16 | 6 | 15 | 0 | 16 | 15 | 10 | 107 | 7.637402 |
| Embiotoca jacksoni | 0 | 14 | 7 | 0 | 13 | 17 | 15 | 12 | 11 | 89 | 6.352605 |
| Paralabrax clathratus | 15 | 8 | 5 | 15 | 9 | 7 | 7 | 7 | 7 | 80 | 5.710207 |
| Embiotoca lateralis | 0 | 0 | 9 | 14 | 11 | 16 | 0 | 6 | 16 | 72 | 5.139186 |
| Oxylebius pictus | 14 | 7 | 0 | 4 | 10 | 13 | 9 | 2 | 7 | 66 | 4.710921 |
| Sebastes atrovirens | 6 | 10 | 10 | 2 | 8 | 9 | 0 | 5 | 12 | 62 | 4.425410 |
| Hypsypops rubicundus | 0 | 9 | 0 | 13 | 7 | 7 | 8 | 9 | 0 | 53 | 3.783012 |
| Rhacochilus vacca | 0 | 5 | 13 | 9 | 0 | 0 | 0 | 4 | 13 | 44 | 3.140614 |
| Hypsurus caryi | 7 | 1 | 6 | 0 | 12 | 0 | 0 | 0 | 9 | 35 | 2.498216 |
| Sebastes caurinus | 9 | 4 | 0 | 6 | 0 | 2 | 0 | 0 | 8 | 29 | 2.069950 |
| Brachyistius frenatus | 10 | 0 | 15 | 0 | 0 | 0 | 3 | 0 | 0 | 28 | 1.998572 |
| Rhinogobiops nicholsii | 0 | 1 | 0 | 0 | 14 | 0 | 0 | 10 | 2 | 27 | 1.927195 |
| Sebastes serranoides | 0 | 3 | 0 | 0 | 0 | 3 | 11 | 8 | 0 | 25 | 1.784440 |
| Hyperprosopon anale | 0 | 0 | 17 | 0 | 0 | 7 | 0 | 0 | 0 | 24 | 1.713062 |
| Girella nigricans | 0 | 12 | 0 | 10 | 0 | 0 | 0 | 1 | 0 | 23 | 1.641684 |
| Others | 32 | 32 | 34 | 42 | 27 | 56 | 59 | 44 | 36 | 362 | 25.838687 |
Cada punto es el “score” de la especie para cada localidad (tabla de arriba).
p1 <- bvi_results %>%
select(-c(BVI, rBVI)) %>%
gather(Sample, Score, -1) %>%
set_colnames(value = c("Spp", "Sample", "Score")) %>%
transform(Spp = reorder(Spp, Score)) %>%
ggplot(aes(x = Spp, y = Score)) +
geom_boxplot(outlier.shape = NULL,
outlier.alpha = 0,
fill = "gray",
alpha = 0.5) +
geom_point(aes(fill = Sample), pch = 21, color = "black", size = 2) +
theme_bw() +
coord_flip() +
labs(x = "Spp", y = "Score") +
scale_fill_brewer(palette = "Paired")
p1
Version animada
plotly::ggplotly(p1)
Podemos visualizar los resultados del IBV (y buscar transiciones en las especies). En este caso, las localidades estan arregladas latitudinalmente, y las especies segun su importancia determinada por el IVB. Las dos figuras son iguales, pero con “smoothing” en la superficie. Sobre todo en el centro, se puede ver como algunas especies desaparecen (al movernos horizontalmente) y otras aparecen. Las de hasta arriba son las constantes, las de hasta abajo son las raras.
bvi_results_plot <- bvi_results %>%
filter(!Spp == "Others") %>%
gather(Location, Score, -c(BVI, rBVI, Spp)) %>%
select(Location, Spp, Score, BVI, rBVI) %>%
left_join(site_data, by = c("Location" = "location")) %>%
transform(Location = reorder(Location, -latitude)) %>%
transform(Spp = reorder(Spp, rBVI))
ggplot(bvi_results_plot, aes(x = Spp, y = Location, fill = Score)) +
geom_raster() +
coord_flip() +
scale_fill_gradientn(colours = colorRamps::matlab.like(20))
ggplot(bvi_results_plot, aes(x = Spp, y = Location, fill = Score)) +
geom_raster(interpolate = T) +
coord_flip() +
scale_fill_gradientn(colours = colorRamps::matlab.like(20))
Ahora hacemos el SIMPER con una aproximacion similar. En este caso, las densidades las transformamos por raiz cuadrada y corremos 999 permutaciones. Al graficar, presentamos solamente las especies que contribuyen al 95% de las disimilitudes. En el eje x presentamos la disimilitud promedio explicada por cada especie.
simper_data <- filter(kelp, transect > 0) %>%
filter(year == 2011,
!location %in% localidades_excluidas) %>%
group_by(location, site, level, transect, genusspecies) %>%
summarize(n = sum(abundance)) %>%
ungroup() %>%
mutate(n = sqrt(n)) %>%
spread(genusspecies, n, fill = 0)
comm <- simper_data %>%
select(-c(location, site, level, transect)) %>%
as.matrix()
sim <- simper(comm = comm,
group = simper_data$location,
permutations = 999,
parallel = 3)
cusums <- tibble(species = sim$ASA_BMA$species) %>%
cbind(map_df(sim, "cusum")) %>%
gather(pair, cumsum, -species)
overall <- map_df(sim, "overall") %>%
gather(pair, overall)
p1 <- tibble(species = sim$ASA_BMA$species) %>%
cbind(map_df(sim, "average")) %>%
gather(pair, average, -species) %>%
left_join(cusums, by = c("species", "pair")) %>%
left_join(overall, by = "pair") %>%
filter(cumsum <= 0.95) %>%
transform(species = reorder(species, average)) %>%
mutate(pair = str_replace(pair, "_", " vs. "),
average = average / overall) %>%
ggplot(aes(x = species, y = average)) +
geom_boxplot(outlier.shape = NULL,
outlier.alpha = 0,
fill = "gray",
alpha = 0.5) +
geom_point(aes(fill = pair),
pch = 21,
color = "black",
size = 2) +
theme_bw() +
coord_flip() +
labs(x = "Spp", y = "Avg % disim")
p1
Version animada
plotly::ggplotly(p1)
Ahora hacemos un ANOSIM incluyendo las localidades de Cedros y San Benito, pero al haber identificado diferencia sen comp, nos quedamos solamente con presencia / ausencia.
anosim_c_cedros <- kelp %>%
filter(year < 2013,
!location %in% localidades_excluidas) %>%
group_by(location, site, level, genusspecies) %>%
summarize(abundance = sum(abundance)) %>%
mutate(abundance = 1) %>%
ungroup() %>%
spread(genusspecies, abundance, fill = 0)
data_c_cedros_groups <- anosim_c_cedros %>%
select(location, site, level) %>%
mutate(loc_site = paste(location, site, sep = "-"))
data_c_cedros_samples <- anosim_c_cedros %>%
select(-c(location, site, level)) %>%
vegdist(method = "bray")
set.seed(43)
mds <- metaMDS(data_c_cedros_samples, trace = F)
stress <- paste("2D Stress =", formatC(mds$grstress, digits = 4, format = "f"))
cbind(data_c_cedros_groups, scores(mds)) %>%
ggplot(aes(x = NMDS1, y = NMDS2, sitio = site)) +
geom_point(size = 4, aes(shape = level, color = location)) +
coord_equal() +
scale_color_brewer(palette = "Paired") +
annotate(geom = "text", x = 0, y = 0.5, label = stress)
perm <- how(nperm = 999)
adonis(data_c_cedros_samples ~ location, data = data_c_cedros_groups)
##
## Call:
## adonis(formula = data_c_cedros_samples ~ location, data = data_c_cedros_groups)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## location 10 2.8044 0.280439 3.8757 0.572 0.001 ***
## Residuals 29 2.0984 0.072359 0.428
## Total 39 4.9028 1.000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
FALTA INCORPORAR ESTOS DATOS, PREGUNTAR A ARTURO POR ELLOS
Hacemos un nMDS comparandoo a nivel de año - localidad - sitio - nivel. En la primer columna se presenta la informacion total. La segunda presenta la transicion de cada localidad / sitio / nivel al conectar los puntos de de 2011 con los de 2013. A y B presentan nMDS con densidades (transformadas por raiz) y C y D presencia / ausencia. Excluimos las muestras de 2012 (Cedros) y los sitios para los que no hay muestras pareadas 2011 - 2013.
anosim_c_cedros <- kelp %>%
filter(location %in% c("ASA", "BMA", "ERE", "ERO", "ISME", "ITSP", "RET", "SMI", "SSI")) %>%
group_by(year, location, site, level, transect, genusspecies) %>%
summarize(abundance = sum(abundance)) %>%
ungroup() %>%
group_by(year, location, site, level, genusspecies) %>%
summarize(abundance = mean(abundance)) %>%
ungroup() %>%
mutate(abundance = sqrt(abundance)) %>%
spread(genusspecies, abundance, fill = 0)
data_c_cedros_groups <- anosim_c_cedros %>%
select(year, location, site, level) %>%
mutate(loc_site = paste(location, site, sep = "-"))
data_c_cedros_samples <- anosim_c_cedros %>%
select(-c(year, location, site, level)) %>%
vegdist(method = "bray")
set.seed(43)
mds <- metaMDS(data_c_cedros_samples, trace = F)
stress <- paste("2D Stress =", formatC(mds$grstress, digits = 4, format = "f"))
p1 <- cbind(data_c_cedros_groups, scores(mds)) %>%
ggplot(aes(x = NMDS1, y = NMDS2)) +
geom_point(aes(shape = level, color = location), size = 4) +
coord_equal() +
scale_color_brewer(palette = "Paired") +
annotate(geom = "text", x = 0, y = 0.5, label = stress)
p2 <- cbind(data_c_cedros_groups, scores(mds)) %>%
ggplot(aes(x = NMDS1, y = NMDS2)) +
geom_line(aes(group = paste(location, site, level))) +
geom_point(aes(shape = paste(year, level), color = location, fill = location), size = 4) +
coord_equal() +
scale_color_brewer(palette = "Paired") +
scale_fill_brewer(palette = "Paired") +
annotate(geom = "text", x = 0, y = 0.5, label = stress) +
scale_shape_manual(values = c(1, 2, 21, 24))
# Reciclamos el codigo y lo hacemos ahora por presencia / ausencia
anosim_c_cedros <- kelp %>%
filter(location %in% c("ASA", "BMA", "ERE", "ERO", "ISME", "ITSP", "RET", "SMI", "SSI")) %>%
group_by(year, location, site, level, genusspecies) %>%
summarize(abundance = sum(abundance)) %>%
ungroup() %>%
mutate(abundance = 1) %>%
spread(genusspecies, abundance, fill = 0)
data_c_cedros_groups <- anosim_c_cedros %>%
select(year, location, site, level) %>%
mutate(loc_site = paste(location, site, sep = "-"))
data_c_cedros_samples <- anosim_c_cedros %>%
select(-c(year, location, site, level)) %>%
vegdist(method = "bray")
set.seed(43)
mds <- metaMDS(data_c_cedros_samples, trace = F)
stress <- paste("2D Stress =", formatC(mds$grstress, digits = 4, format = "f"))
p3 <- cbind(data_c_cedros_groups, scores(mds)) %>%
ggplot(aes(x = NMDS1, y = NMDS2)) +
geom_point(aes(shape = level, color = location), size = 4) +
coord_equal() +
scale_color_brewer(palette = "Paired") +
annotate(geom = "text", x = 0, y = 0.5, label = stress)
p4 <- cbind(data_c_cedros_groups, scores(mds)) %>%
ggplot(aes(x = NMDS1, y = NMDS2)) +
geom_line(aes(group = paste(location, site, level))) +
geom_point(aes(shape = paste(year, level), color = location, fill = location), size = 4) +
coord_equal() +
scale_color_brewer(palette = "Paired") +
scale_fill_brewer(palette = "Paired") +
annotate(geom = "text", x = 0, y = 0.5, label = stress) +
scale_shape_manual(values = c(1, 2, 21, 24))
plot_grid(p1, p2, p3, p4, labels = "AUTO", ncol = 2)
El ANOSIM lo hacemos con densidades promediadas y transformadas por raiz cuadrada.
anosim_c_cedros <- kelp %>%
filter(location %in% c("ASA", "BMA", "ERE", "ERO", "ISME", "ITSP", "RET", "SMI", "SSI"),
transect > 0) %>%
group_by(year, location, site, level, transect, genusspecies) %>%
summarize(abundance = sum(abundance)) %>%
ungroup() %>%
mutate(abundance = sqrt(abundance)) %>%
spread(genusspecies, abundance, fill = 0)
data_c_cedros_groups <- anosim_c_cedros %>%
select(year, location, site, level) %>%
mutate(loc_site = paste(location, site, sep = "-"),
year_loc_site = paste(year, location, site, sep = "-"))
data_c_cedros_samples <- anosim_c_cedros %>%
select(-c(year, location, site, level)) %>%
vegdist(method = "bray")
perm <- how(nperm = 999)
adonis(formula = data_c_cedros_samples ~ year + location, data = data_c_cedros_groups, permutations = perm)
##
## Call:
## adonis(formula = data_c_cedros_samples ~ year + location, data = data_c_cedros_groups, permutations = perm)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## year 1 1.035 1.03547 5.3095 0.02273 0.002 **
## location 8 6.680 0.83497 4.2814 0.14665 0.001 ***
## Residuals 194 37.834 0.19502 0.83062
## Total 203 45.549 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Le falta mucho diseño al mapa, pero es la idea… hacer zoom a las areas de interes
proj <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
data(World)
World <- st_as_sf(World) %>%
st_transform(crs = proj) %>%
filter(iso_a3 == "USA")
coastline <- readRDS(here("raw_data", "spatial", "coastline_mx.rds")) %>%
st_as_sf()
points <- kelp %>%
group_by(year, location, site, latitude, longitude) %>%
count()
major <- ggplot() +
geom_sf(data = World, fill = "gray") +
geom_sf(data = coastline, fill = "gray") +
geom_point(data = points, aes(x = longitude, y = latitude), size = 2, fill = "steelblue", shape = 21)
general <- major +
xlim(c(-120, -110)) +
ylim(c(27, 35))
north <- major +
xlim(c(-117.5, -116)) +
ylim(c(31.5, 32.5))
center <- major +
xlim(c(-116.5, -115)) +
ylim(c(29.5, 31.5))
south <- major +
xlim(c(-116, -114.5)) +
ylim(c(28, 28.5))
left <- plot_grid(general, NULL, labels = c("A", NA), ncol = 1)
right <- plot_grid(north, center, south, labels = c("B", "C", "D"), ncol = 1)
plot_grid(left, right, labels = c("A", NA), ncol = 2)
sessionInfo()
## R version 3.5.0 (2018-04-23)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 16299)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=Spanish_Mexico.1252 LC_CTYPE=Spanish_Mexico.1252
## [3] LC_MONETARY=Spanish_Mexico.1252 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Mexico.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] bindrcpp_0.2.2 forcats_0.3.0 stringr_1.3.0
## [4] dplyr_0.7.4 purrr_0.2.4 readr_1.1.1
## [7] tidyr_0.8.0 tibble_1.4.2 tidyverse_1.2.1
## [10] tmap_1.11-2 sf_0.6-1 cowplot_0.9.2
## [13] ggplot2_2.2.1.9000 magrittr_1.5 here_0.1
## [16] vegan_2.5-1 lattice_0.20-35 permute_0.9-4
##
## loaded via a namespace (and not attached):
## [1] colorspace_1.3-2 deldir_0.1-15 colorRamps_2.3
## [4] class_7.3-14 gdalUtils_2.0.1.14 leaflet_2.0.0
## [7] rgdal_1.2-18 rprojroot_1.3-2 satellite_1.0.1
## [10] base64enc_0.1-3 dichromat_2.0-0 rstudioapi_0.7
## [13] lubridate_1.7.4 xml2_1.2.0 codetools_0.2-15
## [16] splines_3.5.0 R.methodsS3_1.7.1 mnormt_1.5-5
## [19] knitr_1.20 geojsonlint_0.2.0 jsonlite_1.5
## [22] tmaptools_1.2-4 broom_0.4.4 cluster_2.0.7-1
## [25] png_0.1-7 R.oo_1.22.0 rgeos_0.3-26
## [28] shiny_1.0.5 httr_1.3.1 compiler_3.5.0
## [31] backports_1.1.2 mapview_2.3.0 assertthat_0.2.0
## [34] Matrix_1.2-14 lazyeval_0.2.1 cli_1.0.0
## [37] later_0.7.1 htmltools_0.3.6 tools_3.5.0
## [40] coda_0.19-1 gtable_0.2.0 glue_1.2.0
## [43] reshape2_1.4.3 gmodels_2.16.2 V8_1.5
## [46] Rcpp_0.12.16 cellranger_1.1.0 raster_2.6-7
## [49] spdep_0.7-7 gdata_2.18.0 nlme_3.1-137
## [52] udunits2_0.13 iterators_1.0.9 crosstalk_1.0.0
## [55] psych_1.8.3.3 rvest_0.3.2 mime_0.5
## [58] gtools_3.5.0 XML_3.98-1.11 LearnBayes_2.15.1
## [61] MASS_7.3-49 scales_0.5.0.9000 hms_0.4.2
## [64] promises_1.0.1 parallel_3.5.0 expm_0.999-2
## [67] RColorBrewer_1.1-2 yaml_2.1.18 curl_3.2
## [70] geosphere_1.5-7 stringi_1.1.7 jsonvalidate_1.0.0
## [73] highr_0.6 foreach_1.4.4 e1071_1.6-8
## [76] boot_1.3-20 spData_0.2.8.3 rlang_0.2.0.9001
## [79] pkgconfig_2.0.1 bitops_1.0-6 evaluate_0.10.1
## [82] bindr_0.1.1 labeling_0.3 htmlwidgets_1.2
## [85] tidyselect_0.2.4 osmar_1.1-7 plyr_1.8.4
## [88] R6_2.2.2 DBI_0.8 haven_1.1.1
## [91] pillar_1.2.1 foreign_0.8-70 withr_2.1.2
## [94] mgcv_1.8-23 units_0.5-1 RCurl_1.95-4.10
## [97] sp_1.2-7 crayon_1.3.4 modelr_0.1.1
## [100] rmapshaper_0.4.0 KernSmooth_2.23-15 plotly_4.7.1
## [103] rmarkdown_1.9 readxl_1.1.0 grid_3.5.0
## [106] data.table_1.10.4-3 digest_0.6.15 classInt_0.2-3
## [109] webshot_0.5.0 xtable_1.8-2 httpuv_1.4.1
## [112] R.utils_2.6.0 stats4_3.5.0 munsell_0.4.3
## [115] viridisLite_0.3.0